
%%% Weighting the Evidence: A Rank-Dependent Model of Outdoor Recreation, June 2024
%%% This function simulates the ranks for the different harvest outcomes given catch-per-trip and catch-at-length distributions

par_nj=importdata('nbin_param_nj.xlsx'); 
par_ny=importdata('nbin_param_ny.xlsx'); 
par_ct=importdata('nbin_param_ct.xlsx'); 
F_size_NJ=importdata('ProbFlukeSizesNJ2021.txt');
F_size_NJ.data=horzcat(F_size_NJ.data,cumsum(F_size_NJ.data(:,2))); 
F_size_NY=importdata('ProbFlukeSizesNY2021.txt');
F_size_NY.data=horzcat(F_size_NY.data,cumsum(F_size_NY.data(:,2))); 
F_size_CT=importdata('ProbFlukeSizesCT2021.txt');
F_size_CT.data=horzcat(F_size_CT.data,cumsum(F_size_CT.data(:,2))); 

n=1000;
m=100;
kk=0.10; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%NEW JERSEY
b=3;
sizelimit=18;

[r1,c1]=find(F_size_NJ.data==sizelimit);
F_NJ=F_size_NJ.data(r1,3)-kk; 

mu=par_nj(4)+par_nj(1)*0.97+par_nj(2)*8.72+par_nj(3)*8.72^2;
r=1/exp(par_nj(5));
p=r/(r+mu);

%Calculating the probability of keeping zero fish.
sum=0;
for i=1:n
p0=nbinpdf(i,r,p)*F_NJ^i;
sum=sum+p0;
end
p0=sum+nbinpdf(0,r,p);

%Calculating the probablility of keping 1,...,b-1 fish.
P=zeros(1,b-1);
for k=1:b-1
    sum=0;
    for i=k:n
        pk=nbinpdf(i,r,p)*nchoosek(i,k)*F_NJ^(i-k)*(1-F_NJ)^k;
        sum=sum+pk;
    end
    P(k)=sum;
end

%Calculating the probability of keeping the bag limit b.
pb=1-p0-P(1)-P(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%NEW YORK
b=4;
sizelimit=19;

[r1,c1]=find(F_size_NY.data==sizelimit);
F_NY=F_size_NY.data(r1,3)-kk; 

mu=par_ny(4)+par_ny(1)*0.97+par_ny(2)*8.72+par_ny(3)*8.72^2;
r=1/exp(par_ny(5));
p=r/(r+mu);

sum=0;
for i=1:n
p0=nbinpdf(i,r,p)*F_NY^i;
sum=sum+p0;
end
p0=sum+nbinpdf(0,r,p);

P=zeros(1,b-1);
for k=1:b-1
    sum=0;
    for i=k:n
        pk=nbinpdf(i,r,p)*nchoosek(i,k)*F_NY^(i-k)*(1-F_NY)^k;
        sum=sum+pk;
    end
    P(k)=sum;
end

pb=1-p0-P(1)-P(2)-P(3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%CONNECTICUT
b=4; 
sizelimit=19;

[r1,c1]=find(F_size_CT.data==sizelimit);
F_CT=F_size_CT.data(r1,3)-kk; 

mu=par_ct(4)+par_ct(1)*0.97+par_ct(2)*8.72+par_ct(3)*8.72^2;
r=1/exp(par_ct(5));
p=r/(r+mu);

sum=0;
for i=1:n
p0=nbinpdf(i,r,p)*F_CT^i;
sum=sum+p0;
end
p0=sum+nbinpdf(0,r,p);

P=zeros(1,b-1);
for k=1:b-1
    sum=0;
    for i=k:n
        pk=nbinpdf(i,r,p)*nchoosek(i,k)*F_CT^(i-k)*(1-F_CT)^k;
        sum=sum+pk;
    end
    P(k)=sum;
end

pb=1-p0-P(1)-P(2)-P(3);
